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Abstract 

We are interested in bounding probabilities of rare events in the context of com- 
puter experiments. These rare events depend on the output of a physical model with 
random input variables. Since the model is only known through an expensive black 
box function, standard efficient Monte Carlo methods designed for rare events cannot 
be used. We then propose a strategy to deal with this difficulty based on impor- 
tance sampling methods. This proposal relies on Kriging metamodeling and is able to 
achieve sharp upper confidence bounds on the rare event probabilities. The variability 
due to the Kriging metamodeling step is properly taken into account. 
The proposed methodology is applied to a toy example and compared to more stan- 
dard Baycsian bounds. Finally, a challenging real case study is analyzed. It consists 
of finding an upper bound of the probability that the trajectory of an airborne load 
will collide with the aircraft that has released it. 

Keywords: computer experiments, rare events, Kriging, importance sampling, Bayesian 
estimates, risk assessment with fighter aircraft. 
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1 Introduction 

Rare events are a major concern in the reliability of complex systems (Heidelberg, 1995; 
Shahabuddin, 1995). We focus here on rare events depending on computer experiments. 
A computer experiment (Welch et al., 1992; Koehler and Owen, 1996) consists of an 
evaluation of a black box function which describes a physical model, 

y = m, (1.1) 

where y G M and x G where ii^ is a compact subset of M. The code which computes 
/ is expensive since the model is complex. We assume that no more than calls to / 
are possible. The input x are measured with a lack of precision and some variables are 
uncontrollable. Both sources of uncertainties are modeled by a random distribution on E. 
Let X be the random variable. Our goal is to propose an upper bound of the probability: 



(X) < p)) = P(X G Rp) = Px(i?p) , 

where Rp is a subset of E defined by i?p = {x : /(x) < p} and p G M is a given threshold. 
In a crude Monte Carlo scheme, the following estimator of vTp is obtained: 

r(/, Xi;Af,p) 

VTp.iV = , (1.2) 

where r(/, Xi:Ar,p) is defined by 

N 

r(/, Xi.^v, P) = Y1 i]-oo,p[(/(x.)) , (1.3) 
1=1 

and Xi:7v = (Xi, . . . , X^v) is an iV-sample of random variables with the same distribution 
as X. Its expectation and its variance are: 

E(7rp,7v) = P(X G Rp) = TTp , V(7rp,jv) = ^vrp(l - Tip) . 

Since r(/, Xi:Ar, /?) follows a binomial distribution with parameters N and TTp, an exact 
confidence upper bound on vTp: 

F(7rp < 6(r(/, Xi:,v, p),N,a))>l-a, 

is available. 
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Indeed, let T be a random variable which follows a binomial distribution with parameters 
A'^ and p. For any real number a S [0, 1], we can easily show that the upper confidence 
bound b on p: 

¥t{p < b{T,N,a)) > 1 - a 

is such that: 

\ 6 is the solution of equation Ylk=o (fc')^'^(^ ~ b)^~^ = a otherwise 

This upper bound is not in closed form but easily computable. 

In the case where r(/, Xi:7v, /o) = which happens with probability (1 — vTp)^, the (1 — 
a)-confidence interval is [0, 1 — (a)^/^]. As an example, if the realization of r(/, Xi:Ar,p) 
is equal to 0, an upper confidence bound at level 0.9, TTp < 10~^ can be warranted only if 
more than 230,000 calls to / were performed. 

When the purpose is to assess the reliability of a system under the constraint of a limited 
number of calls to /, there is a need for a sharper upper bound on TTp. Several ways to 
improve the precision of estimation and bounding have been proposed in the literature. 

Since Monte Carlo estimation works better for frequent events, the first idea is to 
change the crude scheme in such a manner that the event becomes less rare. It is what 
importance sampling and splitting methods schemes try to achieve. 

For example L'Ecuyer et al. (2007) showed that randomized quasi-Monte Carlo can be used 
jointly with splitting and/or importance sampling. By analysing a rare event as a cascade 
of intermediate less rare events, Del Moral and Garnier (2005) developed a genealogical 
particle system approach to explore the space of inputs E. Cerou and Guyader (2007a, b) 
proposed an adaptive multilevel splitting also based on particle systems. An adaptive 
directional sampling method is presented by Munoz Zuniga et al. (2010) to accelerate the 
Monte Carlo simulation method. These methods can still need too many calls to / and 
the importance distribution is hard to set for an importance sampling method. 

A general approach in computer experiments is to make use of a metamodel which 
is a fast computing function which approximates /. It has to be built on the basis 
of data {/(xi),--- ,/(x„)} which are evaluations of / at points of a well-chosen design 
Dn = {xi, • • • ,x„}. The bet is that these n evaluations will allow the building of more 
accurate bounds on the probability of the target event. 

Kriging is such a metamodeling tool: one can see Santner et al. (2003) and more recently 
Li and Sudjianto (2005); Joseph (2006); Bingham et al. (2006). The function / is seen as 
a realization of a Gaussian process which is a Bayesian prior. 

The related posterior distribution is computed conditionally to the data. It is still a 
Gaussian process whose mean can be used as a prediction of / everywhere on E and the 
variance as a pointwise measure of the confidence one can have in the prediction. 
By using this mean and this variance, Oakley (2004) has developed a sequential method 
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to estimate quantiles and Vazquez and Beet (2009) a sequential method to estimate the 
probabihty of a rare event. Cannamela et al. (2008) have proposed some samphng strate- 
gies based only on a reduced model which is a coarse approximation of / (no information 
about the accuracy of prediction is given), to estimate quantiles. 

In this paper, we also use Kriging metamodeling. Indeed, we assume that / is a 
realization of a Gaussian process F. This Gaussian process is assumed independent of X 
since it models the uncertainty in our knowledge of / while X models a physical uncertainty 
on the input variables. As a consequence, vTp is a realization of the random variable: 



A natural approach consists of focusing on the posterior distribution of Hp which depends 
on the posterior distribution of / given its computed evaluations. A Bayesian estimator 
of Hp can be computed and a credible bound is reachable by simulating realizations of the 
conditional Gaussian process to obtain realizations of Hp. 

We propose another approach which makes use of an importance sampling method the 
importance distribution of which is based on the metamodel. 

The paper is organized as follows: Section 2 describes the posterior distribution of the 
Gaussian process and how to obtain an estimator and a credible bound of Hp. Section 
3 presents our importance sampling method and the stochastic upper bound which is 
provided with a high probability. Finally in Section 4, the two methods are compared on 
a toy example. Different designs of numerical experiments (sequential and non sequential) 
are performed for this comparison. A solution to a real aeronautical case study about the 
risk that the trajectory of an airborne load will collide with the aircraft that has released 
it, is proposed. 

2 Standard Bayesian bounds 

The first step for Kriging metamodeling is to choose a design Dn = {xi, . . . ,x„} of nu- 
merical experiments (one can see Morris and Mitchell (1995); Koehler and Owen (1996) 
and more recently Fang et al. (2006); Mease and Bingham (2006); Dette and Pepelyshev 
(2010)). Let yD„ = (vi = /(xi), . . . ,y„ = /(x„)) be the evaluations of / on 
Let us start from a statistical model consisting of Gaussian processes Fj^ ^ g, the expressions 
of which are given by: for x G E, 



np = E(I]_^,,[(F(X))|F). 



L 





where 
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• hi, . . . ,hL are regression functions, and /3 = . . . , Pl) is a vector of parameters, 

• C is a centered Gaussian process with covariance 

Cov(C(x),C(x')) = a2i^e(x,x'), 

where Kq is a correlation function depending on some parameters (for details 
about kernels, see Koehler and Owen, 1996). 

The maximum likelihood estimates /3, o", of /3, o", are computed on the basis of the 
observations. Then, the Bayesian prior on / is chosen to be F = ^ ^ and the process 

F is assumed independent of X. We denote F^'^ the process F conditionally to i^(xi) = 
yi, . . . , F(x„) = ?/„, in short Yd„ = 

The process F^" is still a Gaussian process (see Santner et al., 2003) with 

• mean: Vx, 

mz5„(x) = H{^)^P + - J) , (2.2) 

• covariance: Vx, x', 

i^^„(x,x') = <72(K^(x,x') - S^^^S^i^^S,,^J , (2.3) 

where 

(SD„D„)i<i,i<n = K0(xi,Xj) and SxD„ = (^^(x, x^))^^.^^ . 

In this approach the conditioning to the data regards the parameters as fixed although 
they are estimated. 

The Bayesian prior distribution P^? on / leads to a Bayesian prior distribution on 
Hp. Our goal is to use the distribution of the posterior process F-^" conditionally to 
the observation of Yo^, to learn about the posterior distribution of lip. It is straight- 
forward to show that Hp given Yo^ = yD„, denoted 11^", has the same distribution as 
E(I]_oo,p[(-F'^" (X))!^-^"). The mean and the variance of 11^" are then given by: 

E(n^") = |^E(I]_^,p[(F^"(x)))Px(dx) = e(^^ ( ^i^^^x X) )) ' ^^-^^ 

where $ is the cumulative distribution function of a centered reduced Gaussian random 
variable, 

V(n^") = / Cov(I]_^,p[(F^"(x)),I]_^,p[(F^"(x'))Px X Px(f^x,dx') . (2.5) 

JExE 
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A numerical Monte Carlo integration can be used to compute the posterior mean and 
variance since they do not need more calls to /. However, the computation time requested 
by a massive Monte Carlo integration, especially for V(n^"), can be very long as can be 
seen in the examples. 

The mean and the variance of 11^" can be used to obtain credible bounds. As a 
consequence of Markov inequality, it holds, for any a G [0, 1], 

IP n^" < " ' \ >\-a. (2.6) 




Likewise, Chebychev inequality gives, for any a G [0, 1], 



n^" <E(n^") + J ^^'"^"^ I > 1-a. (2.7) 



a 



Moreover, the quantiles of 11^" can be estimated through massive simulations of the 
conditional process F^" . These realizations of lead to realizations of H^" from which 
quantiles can be estimated. These quantiles of 11^" are exactly the upper bounds that are 
sought. 

We adapt the algorithm proposed by Oakley (2004) to obtain realizations of H^". 
From a realization , the corresponding realization of 11^" is computed using a massive 
Monte Carlo integration with respect to the distribution of X. Thus, a credible bound on 
n^" is constructed. Given a G (0, 1), a constant a G [0, 1] is found such that: 

P(n^" < a) > 1 - a . 

However, it is not possible to get an exact realization of F^"^ or to sample jointly 
i^^" at all the inputs of the Monte Carlo sample (of X) since this sample is too large. 
Thus, the algorithm relies on a discretization of the process. The same scheme as the 
one followed by Oakley (2004) is used. We choose T points in E: D' = {x'^, . . . , x'j,} 
where the corresponding realizations of F^" are simulated. From the set {y^, . . . ,2/2"} 
joint realizations of {-F^"(x']^), . . . , i<'-^"(x^)}, a realization of is approximated by 
the mean of F^"^'^ which is the process F conditioned to -F(xi) = yi,.. ■ ,F(x„) = y„ 
and F{'k'i) = y'l,. ■ ■ ,-F(x^p) = y'j,. The variance of F^"'^ has to be very small for any 
point in E for the approximation to be valid. Hence, T has to be large enough and the 
points in D' have to fill the space E. When the dimension of the input space is low, 
to propose such a set D' is quite easy. However, it can be burdensome to fill the space 
in high dimension and it can lead to a too large number of needed realizations of F^" 
which are impossible to simulate jointly. The discretization step is a major concern since 
it induces an uncontrollable error on the credible bound on H^" . We propose then in the 
next Section, an alternative approach based on importance sampling which avoids this 
discretization step. 
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3 Metamodel-based importance sampling 



As was explained in Section 1, the major drawback of the crude Monte Carlo scheme is 
the high level of uncertainty when it is used for estimating the probability of a rare event. 
Importance sampling is a way to tackle this problem. The basic idea is to change the 
distribution to make the target event more frequent. We aim at sampling according to 
the importance distribution: 

where Rp C E is to be designed close to Rp = {x £ E : /(x) < p}. Thanks to n calls to 
the metamodel, a set Rp can be chosen as follows: 

Rp = Rp,K. = |x : mD„(x) < /9 + Kv^i^D„(x,x)| , (3.1) 

where k is fixed such that "{x : i^(x) < p} C Rp.n with a good confidence level". In other 
words, if X is such that /(x) < p, we want x to be in Rp^K- We recall that the posterior 
mean rriDni^) approximation of /(x) and Ky^K£)^{:x.,x) has been added to take into 
account the uncertainty of the approximation. 

A set of m points, Zi-m = (Zi, . . . , Zm), is drawn to be an i.i.d. sample following the 
importance distribution. The corresponding importance sampling estimator of Tip is 



m 



T{f, Zi.„) = Yl I]-oo,p[(/(Zfc)) . (3.2) 



m m 

k=l 



The probability Px(-Rp) is computable by a Monte Carlo integration since it does not 
depend on /; yet, m more calls to / are necessary to compute I]_oo,p[(/(Zfc)). This esti- 
mator is only unbiased provided that Rp C Rp. Nevertheless, it is an unbiased estimator 
of Ex(I]-oo.p[(/(^))I^ (^))- Since r(/, Zi:^) follows a binomial distribution 

/ E(I]_^ p[(/{X))I^ (X))\ 

B m, - — 2 ^ for any a g]0; 1[, the following confidence upper bound holds: 



IE(I]-oo,p[(/(X))I^/X))<6(r(/,Zi;^,p),m,a)Px(iip)j >l-a, (3.3) 

by using the bound (1.4). This is an upper bound on vTp only if the estimator (3.2) is 
unbiased i.e. only if Rp C Rp. As is noticed in the decomposition: 

vTp = E(I]_^,,[(/(X))) = E(I]_^,,[(/(X))I^^(X)) + E(I]„^,,[(/(X))(1 - I^^(X))) , 

the second term on the right-hand side which is the opposite of the bias has to be controlled. 
That is why the random variable 
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whose a realization is TTp, is considered. 

Similarly to the previous decomposition, it holds that 

U^" = E(I]_^,,[(F^"(X))I^^(X)|F^") +IE(I]-oo,p[(i^''"(X))(l -I^/X))|F^") . (3.4) 
A bound on E(I]_^_p[(F^"(X))I^^(X)|F^") comes from (3.3). 
Proposition 3.1. For a g]0, l[, it holds that 

P((e(I]_^,,[(F^"(X))I^^(X)|F^") < b Px(/^p)) > 1 - a, (3.5) 

where b stands for b{r{F^", Zi-m, p), fn, a). 
Proof 

Let If be any realization of F^" . 
As in (3.3), we have 

p(e(I]_^,^[((^(X))I^^(X)) < 5(r(</^,Zi:^,p),m,Q)Px(iip)) > 1-a. 
Thus, since this result holds for any realization of F^"- , 

p(e(I]_^,^[(F^"(X))I^^(X)|F^") < 6(r(F^",Zi;^,/5),m,a)Px(iip)) > 

□ 

The next proposition states an upper bound for the second term in (3.4). 
Proposition 3.2. For /3 g]0, 1[, it holds that 

p(e(I]_^,,[(F^"(X))(1-I^^(X))|F^") < ^) > 1-/3, 

where c = E f $ f ^^^dJ^ \ _ j (x))^ . 
Proof 

The mean of E(I]_oo p[(F^"(X))(l — (X))|F^") can be computed in the same fashion 
as the mean of 11^" . It gives 

E (e(I]_,,[(^^^"(X))(1 - I^JX))|F^")) = E (^l^^^ (1 - I^^(X)) 

Then, Markov inequality is applied which completes the proof. □ 
Finally, by gathering the results of Proposition 3.1 and Proposition 3.2, a stochastic upper 
bound is found on 11^" . 
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Proposition 3.3. For a,/3 g]0, 1[ such that a + /3 < 1, it holds that 



n^" <bPx(i?p) + ^j >l-(a + /3), (3.6) 

where b and c have been defined above. 
The proof is obvious. 

If Rp is chosen as proposed in (3.1), the bound c is: 

c = cU)=e('^( ^ -^A,(X) \ ( p-mnA^) 



4 Numerical experiments 
4.1 A toy example 

We study on a toy example, described below, the two bounding strategies: the Bayesian 
strategy (credible bound obtained by simulating realizations of F^" as described in the end 
of Section 2) and the MBIS (metamodel-based importance sampling) strategy (stochastic 
bound given by Proposition 3.3). Since the credible bounds in the Bayesian strategy are 
directly derived from the metamodeling, the Bayesian strategy is the reference and if the 
dimension of E is low, it should perform well. Our aim is to test whether the MBIS 
strategy can achieve such good bounds. Since the choice in the design Dn should directly 
impact the set Rp,K, (3.1) and hence the quality of the importance sampling, different kind 
of designs are considered to compare the strategies. 

The function / : E' = [— 10, 10]^ — )• M+ is assumed to describe a physical model: 

sin(xi) sin(x2 + 2) 
f{xi,X2) = — — + 2 . 

Xi X2 + 2 

The input vector X is supposed to have a uniform distribution on E. The threshold 
is set to p = 0.01 which corresponds to the probability 

Px (/(X) < p) = 4.72 • 10-^. This probability was computed thanks to a massive Monte 
Carlo integration. 

It is assumed that no more than N = 100 calls to / are allowed. For the Bayesian strategy, 
all of the N = 100 available calls to / are used to build the metamodel, while for the MBIS 
strategy (using notations of Section 3) n = 50 and m = 50 are set. The two strategies 
are compared with different design sampling methods. Three design sampling methods 
are used: an LHS-maximin method (Morris and Mitchell, 1995) which is non sequential 
and space filling and two sequential methods. The sequential sampling methods are based 
on a first LHS-maximin design including 80% of the points and the last 20% are added 
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-1 -1 



Figure 1: The function / 

sequentially according to the criterion tIMSE (targeted Integrated Mean Square Error) 
proposed by Picheny et al. (2010) for one method and according to the criterion J SUR 
(Stepwise Uncertainty Reduction) proposed by Beet et al. (2011) for the other method. 
These sequential methods are based on a trade-off between reduction of uncertainty in the 
knowledge of / and the exploration of the space around the critical set Rp. All Kriging 
metamodels are built with an intercept as the regression function and a Gaussian correla- 
tion function is chosen as the correlation function of the Gaussian process C i.e. Vx E E, 
/i(x) = 1 and Vx, x' G E, i('(x, x) = exp (— 0||x — x'|p) are set for the model given by 
equation (2.1). For the Bayesian strategy, a thousand realizations of are computed 
from which the credible bound is obtained. The discretization is done on a grid with 100 
points and to prevent ill-conditioned covariance matrices, if a point of the grid is too close 
to a point of the design it is replaced with a point in E far enough from the points of 
the design and the points of the grid. The numerical integration to compute the realiza- 
tion of from the realization of the process is done with a 10^-sample. In the MBIS 
strategy, the probability Px(-Rp,k) (and also the bound on the bias, given in Proposition 
3.2) was computed by a Monte Carlo integration on a 10^-sample and k = 3 has been set. 

There are sources of variability on the estimators and the bounds due to the design 
sampling methods. Indeed, the designs are computed to be maximin by using a finite 
number of iterations of a simulated annealing algorithm. Moreover, there exist symmetries 
within the class of maximin designs. In the sequential method, the point to be added is 
sought through a stochastic algorithm. Concerning the importance sampling strategy, the 
sampling which gives Zi-m induces variability. In order to test the sensitivity to these 
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sources of variability, each of the two strategies for each of the three design sampling 
methods is repeated one hundred times. 



■3 7 ■ 



Full LHS-maximin 



J SUR 
Design sampling 



tiMSE 



Figure 2: Bayesian 98% credible bound for vTp 
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Minimum 


4.55 


5.50 


4.30 


1*^* quartile 


6.55 


6.20 


6.30 


Mean 


7.78 


7.58 


52.2 


Median 


7.10 


6.30 


6.45 


3'''^ quartile 


7.97 


6.52 


6.87 


Maximum 


35.3 


55.4 


3027 



Table 1: Bayesian 98% credible bound of vTp multiplied by 10^ 

Figure 2 and Table 1 display the results for 98% credible bounds obtained by the 
Bayesian strategy and Figure 3 and Table 2 display the results for 98% stochastic bounds 
obtained by the MBIS strategy. The results are provided according to the design sampling 
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Figure 3: MBIS 98% stochastic bounds of TTp 
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J SUR 


tIMSE 


Minimum 


6.82 


4.98 


5.88 


I''* quartile 


12.58 


5.71 


8.00 


Mean 


16.8 


6.43 


9.67 


Median 


16.4 


6.28 


9.39 


3'''^ quartile 


20.1 


6.81 


11.3 


Maximum 


36.3 


10.4 


16.7 



Table 2: IS 98% stochastic bounds of VTp multiplied by 10^ 



method. For the MBIS stochastic bounds, a = 1% and /? = 1% have been set using the 
notations of Proposition 3.3. If a crude Monte Carlo scheme as presented in Section 1 is 
used here with only N = 100 calls, the estimator is equal to with probability greater than 
0.95 and in this case, the upper confidence bound is 0.038 at level 98%. The two strategies 
bring much sharper bounds on the probability of the rare event. The sequential design 
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sampling method with J SUR criterion leads to the better bounds whatever the strategy. 
The MBIS strategy manages to reach the same sharp bounds as the ones provided by the 
Bayesian strategy in the case where the design is obtained thanks to the J SUR criterion. 
The Bayesian strategy is less sensitive to the choice in the design sampling method. How- 
ever, the Bayesian method suffers from the fact that the quantiles are estimated thanks to 
conditional simulations of the Gaussian process which rely on a discretization of the space. 
Hence, it leads to an approximation in results and stability concern as is noticed in the 
results (see the maximum credible bound obtained with the tIMSE sampling criterion). 
Furthermore, a limited number of iterations is achievable since the simulations are quite 
burdensome. 

The bound provided by Markov inequality (2.6) is not interesting in this example since 
it cannot be less than 0.01. Chebychev inequality has not been used since we were not 
able to determine the posterior variance in a reasonable time. 

As the strategies depend on the Kriging model hypothesis (2.1), a leave-one-out cross 
validation as proposed by Jones et al. (1998) can be performed to check whether this 
hypothesis is sensible. It consists of building n metamodels with posterior mean and 
variance denoted respectively by m^-i and (T^_i, from designs 



where i = 1, . . . ,n. 
Then, the values 



Dn' = {xi,...,Xi_i,Xi+i,...,X„} 

|/(xi) -mj,-.(xi)| 



(4.1) 



are computed. If something like 99.7% of them lie in the interval [—3,3], the Kriging 
hypothesis is not rejected. In our toy example, in all of the tests which were done, all 
these values are in [—2,2]. 

4.2 A real case study: release envelope clearance 
4.2.1 Context 

When releasing an airborne load, a critical issue is the risk that its trajectory could collide 
with the aircraft. The behavior of such a load after release depends on many variables. 
Some are under the control of the crew: mach, altitude, load factor, etc. We call them 
controlled variables and denote their variation domain as C. The others are uncontrolled 
variables: let E be the set of their possible values. The release envelope clearance problem 
consists of exploring the set C to find a subset where the release is safe, whatever the 
uncontrolled variables are. To investigate this problem, we can use a simulator which 
computes the trajectory of the carriage when the values of all the variables are given. 
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Moreover, for xc € C and x £ E, besides the trajectory r(xcx), the program dehvers 
a danger score /(xc,x) to be interpreted as an "algebraic distance": a negative value 
characterizes a collision trajectory. 

To assess the safety of release at a given point of C, we suppose that the values of the 
uncontrolled variables are realizations of a random variable ^ £ E that can be simulated. 
Therefore, for a given value xc G C, and p > the p-collision risk is the probability 

7r^(xc)=P(/(xc,X)</,). 

We do not aim at estimating this risk accurately. 
We would rather classify the points into three categories: according to the position of 
0-risk 7ro(xc) with respect to the two markers 10~^ and 10~^, xc is said to be 

1. totally safe if 7ro(xc) < 10~^, 

2. relatively safe if lO"^ < 7ro(xc) < IQ-^, 

3. unsafe if 7ro(xc) > 10~^. 

In this example, there are 5 controlled and 26 uncontrolled variables, so that C C 
W',E c M'^'\ From budgetary point of view, experts consider that a set of about 400 
representative points of C is enough to cover the domain C consistently. On the other 
hand, the computation of 800, 000 trajectories takes about 4 days which is considered 
reasonable. On the basis of these indications, the maximum amount of available calls to 
the simulator is = 2000 per point. 

4.2.2 Bounding strategy 

As the dimension of the set of uncontrolled variables E is high, the credible bounds ob- 
tained with the Bayesian strategy are impossible to get. The stochastic bounds provided 
by the MBIS strategy are still available. Unfortunately, a sequential sampling method for 
the design of experiments is not achievable since the simulator is too expensive if only 
one point is evaluated per call. Indeed, a fixed part of the cost of a call does not depend 
on the number of points for which the code is run. Although the MBIS strategy with an 
LHS-maximin sampling method is not optimal, it is still efficient. 

We propose this two-step bounding strategy which can be applied for each point of 
the set of representative points (in the set C). Each step uses half of the calls budget: 
m = n = y = 1000. Let xc G C be the current point of interest that we suppose fixed. 
For any x G E, /(x) = /(xc,x) is set. It then corresponds to the notation introduced in 
the first part of the paper. 
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1. At the first stage, a Gaussian process is built as explained in (2), on the basis of 
evaluations /(xi), • • • , /(x„) E M" of / on Dn = {xi, • • • ,x„}. We know that vr^ is 
a realization of the random variable H^" whose mean 



/ p - (X) 



can be computed accurately. 

As stated by (2.6), applying Markov inequality gives, for any a €]0; 1[, 




^" < ^ > 1 - a . 



a 



According to the value of we then take the following decisions: 



p 

-5 



• if E(n^") < ^10^1° which leads by (2.6) to P (nj^^ < > 1 - we 
qualify the current point xc G C as totally safe, 

• if E(n^") > 10~^, we conservatively classify xc as unsafe, 

we use a second stage procedure to refine the risk 

assessment. 

2. A million-sample xi, • • • , xj\,f of X is drawn from which we tune k in such a way that 
m = 1000 of these million elements of E are in Rp^K- The resulting points zi, • • • , 
are an m-sample zi;m of realizations of the random variable Z which follows the 
importance distribution. 



Pz : ^ ^ Px(^|i2p,«) . 

By using m calls to the simulator, r(/, zi-rm p) is computed. Drawn from Proposition 
3.3 with setting a = /3, we obtain the bound 

h{T{f, zi:^, p), m, Q)Px(i?pA) + — > 

a 

which is a decreasing function of a. 

Let us define ao = min{a : 6(r(/, zi:^, p), 'iT', a)Px(-Rp,re) + < 2a}. For such an 
ao, Proposition 3.3 states: 

n^" <6(r(F^",Zi;^,/j),m,ao)Px(i?p) + — ) > 1 - 2qo , 

ao / 

which provides 2ao as a 1 — 2ao confidence upper bound on vTp. 
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4.2.3 Experiments 

Three points of C have been experienced. Of these cases the first one is known to be 
a null 0-risk point, while the third one is very unsafe and the second one is in-between. 
For benchmarking purposes, besides the simulator calls budget required for the estimation 
process described in 4.2.2, a 10, 000 sample of realizations of /(xg, X) has been computed 
for each of the three examples. For each case, we began by estimating a Gaussian process 
on the basis of /-values computed on the points of a 1000-point design Dn = {xi, • • • , x„}. 
This design was obtained by LHS-maximin sampling. Figures 4, 5 and 6 show the predic- 
tive performance of the processes when applied to the benchmark points. These points, 
which appear in red, are sorted according to their process mean values while the blue 
curves mark the predicted 3 standard deviation positions around the means. As appears 
rather clearly, the dispersion of the real values is underestimated by the model: they over- 
flow the blue zone with a frequency (~ 5%) higher than expected (0.27%). The worst case 
is the first one, for which large deviations appear for benchmark points with low values of 
/. In order to obtain bounds from (2.6), we then computed E(n^") using (2.4): 

• In the first case, the massive Monte Carlo procedure leads to a numerically null 
evaluation of E(n^") and, as a consequence, to the classification of the related C 
point as totally safe. 

• In the second example, E(n^") being evaluated at 1.68 10~^, we need to proceed to 
the second step of the bounding strategy to refine the collision probability estimation. 
The obtained confidence upper bound is 1.2 10~^ at confidence level 1 — 1.2 10~^. 
The benchmark data do not show collision 90% confidence upper bound is 
2.3 lO"''. 

• E(n^") = 0.103 in case 3 which is consistent with the 90% confidence interval 
[0.0999; 0.1101], obtained on benchmark data. 

5 Discussion 

In this paper, we have especially focused on bounding the probability of a rare event. 
From our point of view, it seems much more reliable to assess that the probability of a 
feared event (failure of a system, natural disaster, etc.) not exceeding a given level with 
high probability than to estimate the probability of this event happening. Using Kriging 
metamodels to cope with the expensive black box model induces a random interpretation of 
the probability to be estimated. Two strategies were studied in that context including our 
MBIS strategy. On a toy example, the efficiency of the two strategies was shown and it was 
highlighted that a sequential design sampling, when possible, is preferable. Concerning the 
importance sampling strategy, further investigations could be about proposing an optimal 
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Case 1 




Figure 4: Prediction performance case 1 



17 




Figure 5: Prediction performance case 2 
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Figure 6: Prediction performance case 3 
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splitting of the calls to the code used for the metamodel or for the importance sampling. 
Other concerns could be about tuning k to construct the set Rp^K- 

We have dealt with a cross-validation method to assess the Kriging hypothesis. How- 
ever, in the case where the cross-validation leads one to reconsider this hypothesis, a 
solution is to extend the confidence interval on the prediction by tuning by hand the pa- 
rameter cT^ in equation (2.1). In Bayesian words, it can be called using a less informative 
prior distribution on /. 

We have not managed to compute the posterior variance (2.5) by using a massive 
Monte Carlo integration in our examples since it is very small. However, other rare event 
methods can be investigated since the variance no longer depends on /. 
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